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ONR  Project  254921 
Final  report 

I.  Executive  Summary 

Initial  collaborative  observations  and  modeling  of  a  warm  effluent  plume  (Calvert 
Cliffs  Nuclear  Power  Plant)  suggest  that  coherence  length  scales  in  both  models  and 
observations  are  similar  leading  to  confidence  in  the  models  for  AUV  decision 
making. 

Studies  were  initiated  to  determine  if  algorithms  could  be  developed  that  could 
process  current  fields  and  find  regions  of  very  low  current  speed.  Tests  of  the 
algorithm  on  a  two-dimensional  simulated  velocity  field  worked  well  with  true 
stagnation  points  found  and  with  no  ‘false  alarms’. 

II.  Introduction 

This  report  covers  the  initial  project  working  with  colleagues  at  Virginia  Tech  where 
we  made  observations  and  develop  fine  scale  models  in  conjunction  with  their  AUV 
studies.  After  that  we  report  on  progress  working  with  colleagues  at  UDel  on 
stagnation  point  and  center  of  eddy  studies. 

III.  Comparison  of  observations  with  a  model  of  a  warm  water  discharge  into 
Chesapeake  Bay 

Introduction 

The  overall  goal  of  this  project  is  to  develop  and  test  techniques  for  assimilating  data 
from  underwater  and  surface  autonomous  vehicles  in  addition  to  the  usual  sources  of 
Eulerian  and  Lagrangian  systems  into  a  small  scale  coastal  (or  estuarine)  circulation 
model.  The  effort  was  to  focus  on  communication  from  the  various  sensors,  assembly  of 
data  into  gridded  fields  (if  appropriate),  generation  of  nowcasts  and  projection  of 
information  with  forecasts. 

The  project  evolved  to  focus  on  the  comparison  of  the  statistics  of  models  compared  to 
that  of  real  data.  By  doing  this  the  control  systems  may  use  statistical  measures  of  real¬ 
time  data  to  improve  overall  system  control.  The  following  reports  on  this  topic. 

OBJECTIVES 

The  objective  was  to  develop  a  method  to  deterimine  integrative  statistics  such  as 
variance  distributions  and  spatial  correlation  functions  related  to  turbulence  in  the 
thermal  plume.  These  statistics  are  then  used  to  better  control  the  AUVs  in  their  detection 
algorithms,  which  are  different  from  the  gradient  climbing  approaches. 

APPROACH 

The  approach  involves  several  parallel  tasks. 


Assimilative  Modeling  -  The  principle  model  used  was  the  Rutgers/UCLA  Regional 
Ocean  Model  System  (ROMS).  The  NOAA  Coast  Survey  Development  Lab's 
implementation  of  the  QUODDY  finite  element  model  (C3PO)  was  used  to  provide 
velocity  and  sea  surface  height  at  the  model  open  boundaries.  Key  Personnel  were  John 
Klinck,  Mike  Dinniman  (Technical  Staff),  Chester  Grosch,  Andres  Tejada  (Post-Doctoral 
Associate). 

Statistics  and  Spatial  Correlation  Functions  -  Surveys  were  made  in  collaboration 
with  VT  to  gather  data  to  determine  statistics  and  spatial  correlation  functions.  The 
surveys  mainly  were  ADCP  time  series  and  were  spatial  in  nature  to  determine 
spatial  correlation.  Occasional  CTDs  were  used  for  vertical  gradient  determination. 
Key  Personnel  are  Larry  Atkinson,  John  Klinck,  Isaac  Schroeder,  Diego  Naveraz,  Jose 
Blanco  (Post-Doctoral  Associate),  and  David  Salas  (Graduate  Student). 

Model  Configuration 

The  numerical  circulation  model  used  was  the  Rutgers/UCLA  Regional  Ocean  Model 
System  (ROMS)  (Shchepetkin  and  McWilliams,  2005;  Haidvogel  et  al,  2008).  ROMS 
is  a  free-surface,  finite-difference,  hydrostatic  primitive  equation  ocean  circulation 
model  with  a  terrain  following  vertical  coordinate.  A  conservative  parabolic-spline 
discretization  was  implemented  in  the  vertical  in  order  to  reduce  the  pressure 
gradient  error  that  can  often  be  a  problem  in  terrain-following  coordinate  models 
(Shchepetkin  and  McWilliams,  2003).  Momentum  advection  was  computed  with  a 
3rd-order  upstream-biased  advection  scheme  (Shchepetkin  and  McWilliams,  1998). 
Tracers  (temperature  and  salinity  in  this  case)  were  advected  with  a  monotonized, 
fourth-order,  centered  scheme  (Shchepetkin  and  McWilliams,  2003).  Vertical 
mixing  was  simulated  with  the  k  profile  parameterization  (KPP)  vertical  turbulence 
closure  scheme  (Large  et  al.,  1994)  including  surface  and  bottom  (Durski  et  al., 
2004)  boundary  layer  mixing.  No  explicit  horizontal  mixing  of  momentum  or 
tracers  were  used  Quadratic  bottom  stress  with  a  coefficient  of  1.0x10  3  (non- 
dimensional)  was  applied. 

The  model  domain  (Figure  1)  covered  a  3  km  by  4  km  area  of  the  west-central 
Chesapeake  Bay  centered  on  the  cooling  water  discharge  channel  from  the  Calvert 
Cliffs  Nuclear  Power  Plant  located  near  Lusby  in  Calvert  County,  Maryland.  The 
model  grid  was  oriented  to  be  along  the  axis  of  the  discharge  channel.  The 
horizontal  resolution  was  a  constant  15  m  and  the  model  had  20  vertical  layers  that 
were  uniformly  spread  throughout  the  water  column  if  the  bathymetry  was  less 
than  5  m  and  then  slightly  concentrated  towards  the  surface  and  bottom  for  deeper 
water.  One  experiment  was  also  performed  with  a  10  m  horizontal  resolution  grid. 
The  bathymetry  was  from  the  National  Geophysical  Data  Center  3-arc-second 
Coastal  Relief  Model  with  some  modifications  made  for  the  intake  and  discharge 
channels.  In  the  model  domain,  the  bathymetry  ranged  from  2  to  14  m  leading  to  a 
maximum  depth  of  the  model  surface  layer  of  0.46  m. 

Model  initial  conditions  were  based  on  observations  in  the  area  in  mid-August 
(Lacy,  1979)  with  a  constant  temperature  of  26.7  C  and  a  vertical  salinity  profile 
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that  was  a  uniform  9.8  psu  above  5.8  m  and  then  increased  linearly  below  this  at  a 
rate  of  0.45  psu/m.  Wind  stress  was  applied  based  on  a  constant  wind  of  4.1  m/s  (8 
knots)  from  due  south  (roughly  equivalent  to  the  observed  wind  during  the  field 
observations).  Surface  heat  and  fresh  water  fluxes  were  set  to  zero. 

Tidal  sea  level  and  depth  average  velocity  were  imposed  on  the  model  open 
boundaries  using  the  five  most  significant  harmonics  (M2,  01,  N2,  Kl,  S2)  from  the 
NOAA  Coast  Survey  Development  Lab's  Chesapeake  3-D  Physical  Oceanographic 
Model  (C3P0  -  http://nauticalcharts.noaa.gov/csdl/op/c3po.html)  simulation  of  the 
bay  using  the  QUODDY  (Lynch  et  al.,  1996)  finite-element  model.  Note  that  other 
than  the  discharge  jet,  the  observed  predominant  currents  in  the  area  under  typical 
wind  conditions  are  the  tidal  currents  (Lacy,  1979). 

A  two-dimensional  radiation  scheme  was  used  at  the  open  boundaries  for  the 
baroclinic  velocity  components  and  the  tracers  with  adaptive  nudging  also  used  for 
the  tracers  (Marchesiello  et  al.,  2001).The  time  scale  of  the  boundary  nudging  was 
set  to  1  day  and  the  tracers  were  relaxed  to  the  initial  values. 

The  power  plant  has  two  nuclear  units  and  has  a  once-through  cooling  system 
which  draws  in  water  from  the  bay  below  a  curtain  wall  at  the  plant  intake.  At  full 
plant  load,  the  cooling  water  is  heated  between  5  and  6  C  (USAEC,  1973;  Lacy, 

1979).  This  water  is  discharged  through  four  3.8m  by  3.8m  concrete  conduits  that 
rest  on  the  bay  bottom  about  250  m  from  the  shore.  The  tops  of  the  discharge 
conduits  are  about  1.8  m  below  the  water  surface  and  the  velocity  through  each 
conduit  is  about  2.7  m/s  (USAEC  1973;  Schreiner  et  al.,  2002),  leading  to  a  total 
volume  flux  of  156  m3/s  from  the  discharge  channel.  The  discharge  channel  was 
implemented  in  the  model  as  two  point  sources  of  volume  flux  (78  m3/s  each)  that 
are  each  one  grid  cell  wide  and  in  the  vertical  have  flux  into  the  water  from  1.72  m 
below  the  mean  surface  to  the  bottom  (4.30  m).  The  temperature  of  the  discharge 
water  is  set  to  a  constant  32.0  C  (5.3  C  above  ambient)  and  the  salinity  is  set  to  9.8 
psu.  The  velocity  of  the  discharge  jet  at  the  conduit  mouth  is  smaller  than 
observations  (2.0  m/s  vs.  2.7  m/s),  but  it  was  thought  to  be  more  important  to  have 
the  top  location  of  the  plume  with  respect  to  sea  level,  the  total  discharge  volume 
flux  and  the  total  discharge  heat  flux  be  accurate  (having  all  four  quantities  match 
would  require  a  drastic  modification  of  the  bathymetry  in  the  area).The  model 
simulation  was  run  for  36  hours  and  model  fields  were  saved  every  20  minutes. 

Figure  1  shows  the  results  of  a  run  using  similar  tide  and  wind  parameters  as  occurred 
during  the  observation  period. 


Results  of  Model  Runs 


The  Rutgers/UCLA  Regional  Ocean  Model  System  (ROMS)  is  being  used  to  create  a 
model  for  a  region  (3.9  km  by  3.0  km)  in  the  Calvert  Cliffs  area.  This  will  be  the  grid  for 
the  future  assimilative  modeling,  but  is  being  used  now  for  some  preliminary  studies  of 
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the  high  temperature  plume  of  cooling  water  discharged  from  the  Calvert  Cliffs  plant. 

The  model  grid  has  high  horizontal  resolution  (15  m)  and  has  20  vertical  levels.  The 
bathymetry  for  the  region  was  taken  from  National  Geophysical  Data  Center  Coastal 
Relief  Model.  Some  modifications  were  made  by  hand  to  account  for  the  intake  and 
discharge  channels.  The  model  is  forced  winds  at  the  surface  and  by  velocities  and 
surface  heights  at  the  lateral  open  boundaries  representing  the  tidal  flow  in  Chesapeake 
Bay.  There  are  currently  no  surface  heat  or  salt  fluxes.  Figure  1  shows  the  results  of  a  run 
using  similar  tide  and  wind  parameters  as  occurred  during  the  observation  period. 
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SST  Delta  from  Ambient 


Model  day:  0.97222222 


-76.46  -76.45  -76.44  -76.43  -76.42  -76.41 


0.0  1.E5  2.50  3.75  5.00 

deg  C 


Figure  1 .  ROMS  model  simulation  of  the  thermal  plume  emanating  from  the  power  plant 
at  tide  stage  similar  to  when  data  was  taken  (Figure  2).  No  contours  below  1 ,0°C  above 
ambient  are  shown. 
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Field  survey  completed  in  2005  -  A  one-day  survey  of  the  thermal  plume  was 
accomplished  in  early  September.  A  towed  ADCP/CT  system  and  occasional  CTD 
profiles  were  used  to  characterize  the  discharge  area  and  the  plume  both  near  and  far 
field.  Additionally  Dr.  Stilwell  and  students  came  on  board  to  observe  our  techniques  for 
acquiring  data  and  real  time  processing.  This  was  a  very  valuable  experience  for  all 
involved.  Figure  2  shows  one  of  the  surface  SST  fields  derived  from  our  survey. 
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Figure  2.  Surface  temperature  derived  from  towed  ADCP/CT  system.  The  source  of  the 
thermal  plume  is  the  darkest  red.  The  tide  was  flooding  sending  the  plume  northward. 

Realtime  data  needs  assessed  -  Although  we  originally  proposed  two  ADCP  moorings  to 
provide  velocity  data  for  the  model  further  testing  and  model  runs  suggest  we  would  not 
need  the  mooring  data. 
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Statistical  comparison  of  model  and  real  data 


A  key  test  of  the  model  is  a  comparison  of  it  with  the  ADCP  data.  Figure  3  shows  data 
from  the  observations  and  Figure  4  shows  data  from  the  model  run. 


Autocorrelation  of  Surface  Temperature  Anomalies 


Mean  of  the  5  Transacts 


Figure  3.  Five  realizations  of  the  autocorrelation  function.  These  were  calculated  from  the 
temperature  data  obtained  during  five  passes  across  the  outflow.  The  spread  in  the  results 
is  to  be  expected  because  of  the  turbulent  structure  of  the  outflow  jet. 
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Autocorrelation  Autocorrelation 
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Figure  4.  Autocorrelation  from  the  9  transects  in  the  model. 


Figure  4  shows  nine  realizations  of  the  autocorrelation  function  obtained  from  the  model. 
In  both  figures  the  dashed  (upper  panel)  and  the  solid  (lower  panel)  are  the  average  of 
these  realizations.  There  is  reasonable  agreement  between  the  results  obtained  from  the 
data  and  those  from  the  model.  However,  the  first  zero  crossing  of  the  data 
autocorrelation  is  at  approximately  60m  while  that  from  the  model  is  at  about  100m. 

This  indicates  that  there  is  probably  somewhat  too  much  damping  in  the  model. 
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IV.  Stagnation  Point  and  Center  of  Eddy  Studies 


The  original  proposal  called  for  emphasis  on  field  work  providing  data  for  models 
that  would  interface  with  the  research  of  Dr.  Stilwell's  group  at  VT. 


We  did  one  field  experiment  that  was  successful  although  not  in  ways  expected.  We 
implemented  a  very  fine  scale  ROMS  model  of  the  Calvert  Cliffs  nuclear  plant  plume 
that  was  then  compared  to  our  ADCP  data.  We  found  surprisingly  good 
comparisons. 

Further  discussions  with  Dr.  Stilwell  and  new  discussions  with  Drs.  Kirwan  and 
Lipphardt  at  University  of  Delaware  led  us  to  take  an  approach  that  interfaces  with 
both  Dr.  Stilwell  and  Drs.  Kirwan  and  Lipphardt. 

For  the  final  year  of  our  research  on  this  project  we  continued  our  collaboration 
with  Dr.  Stilwell  and  students  and  add  collaboration  with  Drs.  Kirwan  and 
Lipphardt.  This  will  make  the  most  effective  use  of  our  research  skills  and  may  add 
useful  research  to  the  University  of  Delaware  effort. 

We  proposed  to  focus  on  algorithms  which  can  process  the  current  field  with  little 
or  no  human  intervention  and  find  regions  of  very  low  current  speed.  It  appears 
that  the  regions  of  very  low  current  speed,  away  from  boundaries,  must  fall  into  one 
of  two  classes  (1)  near  stagnation  points,  and  (2]  near  the  centers  of  eddies.  Dr. 
Grosch  has  been  working  on  algorithms  which  can  detect  one  or  the  other.  The  next 
step  is  to  test  these  algorithms  on  flows  for  which,  in  a  sense,  the  answer  is  known. 

In  addition,  there  may  be  better  algorithms.  The  results  of  his  collaboration  with 
Drs.  Kirwan  and  Lipphardt  are  as  follows: 

In  certain  circumstances  it  is  desirable  for  an  AVU  to  remain  within  a  given 
volume  with  minimum  power  expenditure.  One  obvious  stategy  is  to  find  regions 
within  the  volume  where  the  velocity  is  very  small,  that  is  near  stagnation  points.  It 
is  necessary  to  develop  one  or  more  algorithms  which  can  be  used  on  velocity  fields 
obtained  from  a  model  and/or  observations.  The  algorithm  should  be  robust  and 
easy  to  implement  without  human  intervention. 

A  possible  algorithm  is  being  implemented  and  tested.  This  algorithm  can  be 
described  as  follows.  Choose  a  small  magnitude,  e ,  not  exactly  zero  to  define  a  near 
stagnation  point.  Assuming  that  the  velocity  data  is  given  on  some,  not  necessarily 
uniform  grid,  step  across  the  grid  computing  the  velocity  magnitude  at  each  point 
and  comparing  it  to  e .  If  the  velocity  magnitude  is  less  than  e  that  grid  point  is  a 
near  stagnation  point. 

In  the  more  likely  situation  that  a  near  stagnation  point  is  not  located  at  a 
grid  point  an  extrapolation  procedure  into  space  near  a  grid  point  can  be  used.  Let 
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«,,/  =  1,2,3  denote  the  velocity  field.  The  velocity  at  a  gridpoint  (x,'0),X2°',xj0))  is 
(m,(0),M20),M30)) .  An  approximation  to  the  velocity  components  in  the  vicinity  of 

(xfVf  ,x<°>)  is 


with  summation  over  j  and  Sfj  the  rate  of  strain  tensor  at  (xf^.xf^xf^ . 

du  du 

If  a  near  stagnation  point  lies  in  the  vicinity  of  (x,'0’,X2°',xj0)) ,  then  its  location  can  be 
found  from 

One  can  check  that  this  is  a  valid  near  stagnation  point  by  testing  whether  or  not 
the  location  determined  is  within  a  sphere  of  radius  R  of  (x^x^xf')  with  R  a 

fraction  of  the  distance  to  the  nearest  neighbor  of  (x^'.x^xf1) . 

This  algorithm  was  coded  and  tested,  first  on  a  two-dimensional  analytic 
velocity  field  which  has  a  stagnation  point,  next  on  a  two-dimensional  simulated 
velocity  field.  The  two-dimensional  simulated  velocity  field  was  a  slice  thru  a  Direct 
Numerical  Simulation  of  decaying  isotropic  turbulence  at  a  Reynolds  number  based 
on  the  Taylor  microscale  of  about  50.  The  algoritm  worked  very  well,  true 
stagnation  points  were  found  and  there  were  no  "false  alarms". 

As  the  project  has  progressed  the  areas  of  fruitful  research  has  changed  also.  One  of  those 
new  areas  has  been  adapting  models  to  optimize  the  locational  control  of  an  AUV. 

The  good  comparison  between  the  model  and  data  suggests  that  we  are  on  the  right  track 
to  improving  the  control  algorithms  via  statistical  analysis  of  the  real  time  data  acquired 
by  the  AUV. 
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